/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined(_Cauchy_p_system_FD_)
#define _Cauchy_p_system_FD_
#define NUMBER_OF_PARTICLES 10
#include <src/math/Operators/all.h>
#include <src/algorithm/Evolution_Algorithm/Crank_Nicolson.h>
#include <src/algorithm/Evolution_Algorithm/Escheme.h>
#include <src/math/CauchyProblems/CauchyProblem_traits.h>
#include <src/math/Operators/Tadmor.h>
#include <src/math/functor/GaussianFunctor.h>



template <class data>
struct EScheme_p_system
	: public EScheme< Wave_visc_t<data> >
{
	typedef typename EScheme_p_system			type;
	typedef  EScheme							Base;
	typedef typename data::value_type			value_type;
	typedef typename data::base_value_type		base_value_type;
	template <class data>
	struct sigma : public integrable<data,data>
	{
		struct int_g : public Base
		{
			virtual result_type operator()(argument_type from) const {
				return pow(from,4)/4.-pow(from,2)/2.; // hyperbolic / elliptic equation
//				return pow(from,2)/2.; // wave equation
			};
		};
		virtual result_type operator()(argument_type from) const {
			return pow(from,3)-from; // hyperbolic / elliptic equation
//			return from;   // wave equation
		};
		virtual smart_ptr<Base> get_integral() const {
			return smart_ptr<Base>(new int_g);
		};
	};
	type() : sigma_( smart_ptr< sigma<base_value_type>> (new sigma<base_value_type> )){
		u_.resize(2,0,0);
	} 
	virtual void vstar(indice_type n) 
	{
		vstar_[n][0] = sigma_(u_[n][0],u_[n+1][0]);
		vstar_[n][1] = (u_[n][1]+u_[n+1][1])/2.;
	};
	virtual void gstar(indice_type n) 
	{
		gstar_[n] = Op_(0,vstar_[n]);
	};
	virtual OS_double Estar(indice_type n) const {
		typedef sigma<value_type> f_type;
		typedef typename smart_ptr<f_type::Base> f_int_type;
		sigma<value_type> f;
		f_int_type f_int_ = f.get_integral();
		return jmmath::norm1((*f_int_)(u_[n][0] )) + jmmath::norm22(u_[n][1])/2.;
	};
	virtual OS_double CFL(indice_type n) {
		OS_double breakpoint1 = jmmath::norm_infty(u_[n]);
		sigma<value_type> f;
		OS_double breakpoint2 = jmmath::norm_infty(f(u_[n][0]));
		OS_double breakpoint3 = std::sqrt(breakpoint2/breakpoint1);

		return .5/ get_state_variable()->size(1)/breakpoint3;

	};
	virtual bool criteria(indice_type n) {
		OS_double breakpoint1 = jmmath::norm22(gstar_[n] - Op_S_k1);
		OS_double breakpoint2 = jmmath::norm22(gstar_[n]);
		return breakpoint1/breakpoint2 < 1e-4;
/*		OS_double breakpoint1 = Estar(n+1); // these lines test the Energy conservation without viscosity (Entropy Conservation)
		OS_double breakpoint2 = Estar(n);
		return (breakpoint1-breakpoint2)/(delta_t_*breakpoint2) < 1e-4;*/
	};
	typedef Tadmor<value_type> tadmor;
	tadmor sigma_;
};

//! Third order scheme for the p system equation \f$\partial_t u = \partial_x v\f$, \f$\partial_t v = \partial_x u\f$.
/*!
    This class initialize the Cauchy problem for the one dimensional p-system equation
    \f[
        \partial_t w = \partial_x z; \partial_t z = \partial_x \sigma(w)
    \f];
    with \f$x \in [0,1]\f$. This system can be written under the scalar form \f$\partial_t^2 w - \partial_x^2 \sigma(w) =0\f$.
    The previous equation owns many weak solutions. To select one ,
    we will consider entropic solutions, i.e. solutions that verifies weakly \f$\partial_t E(u) + \partial_x G(u) \le 0\f$,
    with \f$u:=(w,z)\f$ and
    \f[
        E(u) = \int_{[0,1]} e(w)+ \frac{z^2}{2} := \int_{[0,1]} \left( \int_0^w \sigma(s)ds+ \frac{z^2}{2}
        \right);G(u) = \sigma(w)z
    \f]
    In this setting, the entropy variable is \f$v(u) :=(v_1,v_2)(u) = \nabla_u E(u) =
    (\sigma(w),z)\f$, and the equation may be written in the entropy
    variable as
    \f[
        \partial_t u = g(v) := \left( \partial_x v_2(u), \partial_x v_1(u) \right)
    \f]
    We illustrate the behavior of the finite difference scheme
    \f[ \label{p_system_e}
        \delta_t w^n = \delta_{x} v_2^{*,n+1/2} ; \delta_t z^n = \delta_{\bar{x}} v_1^{*,n+1/2} = 0
    \f]
    where the discrete entropy variable \f$v^{*,n+1/2} \f$
    \f$:= \left( v_{1}^{*,n+1/2},v_{2}^{*,n+1/2} \right)\f$
     \f$:= \{v_{1,i}^{*,n+1/2},v_{2,i}^{*,n+1/2}\}_i\f$ is computed
    accordingly to the (Tadmor) formula
    \f[ \label{p-scheme}
         v_{1}^{*,n+1/2} := \sigma^{*,n+1/2} := \{\frac{\int_{w_i^n}^{w_i^{n+1}}
         \sigma(y)dy}{w_i^{n+1}-w_i^{n}}\}_i; v_{2}^{*,n+1/2}= z^{n+1/2}
    \f]
    This scheme is a second order one, as can be seen from its
    equivalent form
    \f[
        \delta_{t\bar{t}} w^n = \delta_{x} \delta_{\bar{t}}z^{n+1/2}
        = \delta_{x} \delta_{t} \frac{z^{n}+ z^{n-1}}{2} =
        \delta_{x\bar{x}} \frac{\sigma^{*,n+1/2} + \sigma^{*,n-1/2}}{2}
    \f]
    and is conservative for the discrete entropy
    \f[
        {\EE}^n = \sum_i e(w_i^n)+\frac{1}{2} |z_i^n|^2
    \f]
    Let us check that this entropy conservation holds :
    \f[
        \delta_t \frac{1}{2} |z^n|^2 = z^{n+1/2} \cdot \delta_t z^{n} = z^{n+1/2} \cdot
        \delta_{\bar{x}} \sigma^{*,n+1/2}
    \f]
    and
    \f[
        \delta_t e(w^n) = \sigma^{*,n+1/2} \cdot \delta_t w^n =
        \sigma^{*,n+1/2} \cdot \delta_{x} z^{n+1/2}
    \f]
    A discrete integration by parts ends the proof of Entropy conservation.

    Note that this entropy functional may provide a bound if \f$e(u)\f$ is bounded from above, providing a stability result for the scheme in this case.


    Let us turn out to the definition of the CFL condition. The scheme (\ref{p-scheme}) is solved by a fixed point algorithm defined in the class Escheme. This class initiate \f$z^{0,n+1}:=z^n;w^{0,n+1}:=w^n\f$, then loop for
    \f$k\ge1\f$, defining first
    \f[
        z^{k-1,n+1/2} := \frac{z^{k-1,n+1}+z^n}{2};\sigma^{*,k-1,n+1/2} = \{\frac{\int_{w_i^n}^{w_i^{k-1,n+1}}
         \sigma(y)dy}{w_i^{k-1,n+1}-w_i^{n}}\}_i
    \f]
    and compute the \f$k\f$-th iterate \f$(w^{k,n+1},z^{k,n+1})\f$ accordingly to
    \f[
        \delta_t w^{k,n} := \frac{w^{k,n+1}-w^n}{\tau^n} = \delta_x z^{k-1,n+1/2};\delta_t z^{k,n} := \frac{z^{k,n+1}-z^n}{\tau^n} = \delta_{\bar{x}}
        \sigma^{*,k-1,n+1/2}
    \f]
    We may take the scalar product of the right hand side of the previous equation by \f$z^{k,n+1/2}\f$, obtaining with a discrete part integration
    \f[
        \frac{|z^{k,n+1}|^2-|z^{n}|^2}{\tau^n}=z^{k,n+1/2}\cdot \delta_{\bar{x}}\sigma^{*,k-1,n+1/2}=- \delta_{x}z^{k,n+1/2}\cdot \sigma^{*,k-1,n+1/2}
    \f]
    Using the scheme, the previous expression is equal to
    \f[
       -\delta_t w^{k+1,n} \cdot \sigma^{*,k-1,n+1/2}= -\frac{e(w^{k-1,n+1})-e(w^n)}{\tau^n}-\frac{w^{k+1,n+1}-w^{k-1,n+1}}{\tau^n}\cdot \sigma^{*,k-1,n+1/2}
    \f]
    Defining the discrete iterated entropy \f${\EE}^{k,n+1} = \sum_i e(w_i^{k-1,n+1})+\frac{1}{2} |z_i^{k,n+1}|^2\f$, we get
    \f[
       \frac{{\EE}^{k,n+1}-{\EE}^{n}}{\tau^n} = -\frac{w^{k+1,n+1}-w^{k-1,n+1}}{\tau^n}\cdot \sigma^{*,k-1,n+1/2}
    \f]
    that we can write as
     \f[
       \frac{{\EE}^{k,n+1}-{\EE}^{n}}{\tau^n} = -\{\frac{w_i^{k+1,n}-w_i^{k-1,n}}{w_i^{k-1,n}-w_i^{,n}}\}_i\cdot \frac{e(w^{k-1,n+1})-e(w^{n})}{\tau^n}
    \f]
    The first iterate \f$k=1\f$ is
    \f[
       \frac{{\EE}^{1,n+1}-{\EE}^{n}}{\tau^n} = -\frac{w^{2,n}-w^{n}}{\tau^n}\cdot \sigma^{n}=-\delta_x z^n \cdot \sigma^{n}
    \f]
    Entropy conservative schemes produces spurious oscillation at shock formation. To remove them,
    we consider the following (entropy) viscous approximation of p-system
    \f[
        \partial_t w = \partial_x z; \partial_t z = \partial_x \sigma(w)+\epsilon \partial_{xx}v
    \f]
    where the \f$\epsilon\f$ is a small parameter and \f$v\f$ denotes the entropy variable. This leads to consider the viscous approximation of the scheme (\ref{p_system_e})
    \f[
        \delta_t w^n = \delta_{x} v_2^{*,n+1/2} ; \delta_t z^n = \delta_{\bar{x}} v_1^{*,n+1/2} + \frac{1}{N}\delta_{\bar{x}x} v_1^{*,n+1/2}= 0
    \f]
    for which the following discrete entropy dissipation holds
     \f[
       \delta_t \sum_{i=1,..,N} e(w_i^n)+\frac{1}{2} |z_i^n|^2 = \frac{{\EE}^{n-1}-{\EE}^n}{\tau^n} \le 0
    \f]

    The boundary conditions are given by <data> : either Dirichlet or Periodic ones.
    The initial data are \f$u^0=(\delta_x w^0,z^0)\f$, where \f$z^0 \equiv 0\f$ and \f$w^0\f$ is a normal distribution centered in
    0.5. We considered a system of elliptic / hyperbolic type \f$\sigma(w) = w^3-w\f$.

*/
template <class data, class Interface = Graphical_Interface>
		class Cauchy_Entropic_p_system : public CauchyProblem_traits< EScheme_p_system < data >, Interface >
		{
		public :
			virtual ~Cauchy_Entropic_p_system(){}
			Cauchy_Entropic_p_system(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Entropic_p_system);
			};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"Cauchy_Entropic_p_system : only dim =1");
			};
			OS_size get_dim() {return 1;};
			virtual void initial_conditions(OS_size nb_part)
			{
				data result(2,nb_part);
				result[0] = Gaussian_Functor<data::value_type>()(nb_part);
//				result[0] = (result[0]- (result[0] <<1))*nb_part;
				result[0] /= .25*jmmath::norm_infty(result[0]);
//				result[1].assign(0.);
				result[1] = result[0];
				result[1] <<=1;
				set_initial_conditions(result);
			};
			virtual std::string get_name() {return "p-system : " + algorithm_type::get_name()
													+time_operator_type::Boundary_type::get_name();
			};

		};

#undef NUMBER_OF_PARTICLES

#endif

